!************************************************************
!------------------------------------------------------------
! $Id: x2y2model.shar,v 1.2 2011/05/14 21:34:35 daw Exp $
! Description: does MD by velocity verlet
! Author:  Murray S. Daw <daw@clemson.edu>
! Created: 6 Nov 2007
!-------------------------------------------------------------
! Philip Ladd contributed some lines in 2006
! as a mod to original by msd (using RK method)
! Velocity verlet is sympletic and time-reversible
!	q(n+1) = q(n) + p(n)dt + 0.5f(q(n))dt^2
!	p(n+1) = p(n) + 0.5(f(q(n)) + f(q(n+1))dt
!--------------------------------------------------------------
! Modify MD part with new potential
! Author: Yang Gao
! data: 09/23/2012
!--------------------------------------------------------------
module md_vv

! Load modules
  use prec
  IMPLICIT NONE

  real(r8) :: lambda,lambda1              ! (anharmonicity)
  real(r8), dimension(2) :: fs, p,q   ! forces
  real(r8), dimension(2) :: fsnext  ! forces (next step)
  real(r8) :: dt, mr ! timestep (set in setupmd1dvv)

  private :: fs, fsnext, p, q, dt, lambda, lambda1

contains

  subroutine setupmdvv( lambdain, lambda1in, dtin )

    real(r8), intent(in) :: lambdain
    real(r8), intent(in) :: lambda1in
    real(r8), intent(in) :: dtin

! store lambda
    lambda = lambdain
    lambda1 = lambda1in
    dt = dtin

    return
  end subroutine setupmdvv

! ps is momentum NOT velocity
  subroutine mdvv( qs, ps, nsteps )

    real(r8), dimension(2), intent(inout) :: qs
    real(r8), dimension(2), intent(inout) :: ps
    integer, intent(in) :: nsteps
    real(r8):: mom0, momf,e0, ef
    real(r8), dimension(2) ::qsn, psn, pupdate
    integer :: n, INFO
    real(r8):: EPSFCN, TOL, WA(1), LWA


! Declare local parameters
    integer :: istep        ! index for summations

! calculate forces in cartesian coordinate

    fs(1) = -(qs(1)+lambda*(qs(1)**3+qs(1)*qs(2)**2)+ &
             lambda1*(qs(1)**2-qs(2)**2))

    fs(2) = -(qs(2)+lambda*(qs(2)**3+qs(2)*qs(1)**2)- &
             2*lambda1*qs(1)*qs(2))
    
	e0 = (qs(1)**2 + qs(2)**2 + ps(1)**2 + ps(2)**2)/2. + &
	      0.25*lambda*(qs(1)**2 + qs(2)**2)**2 + &
              1./3.*lambda1*(qs(1)**3-3*qs(1)*qs(2)**2)

	mom0 = qs(1)*ps(2)-qs(2)*ps(1)

! Numerically integrate using Velocity Verlet
    do istep = 1,nsteps


	qsn(1) = qs(1) + ps(1) * dt + 0.5*fs(1) * dt**2
	qsn(2) = qs(2) + (ps(2)) * dt + 0.5*(fs(2)) * dt**2

		qs = qsn
	fsnext(1) = -(qs(1)+lambda*(qs(1)**3+qs(1)*qs(2)**2)+ &
                lambda1*(qs(1)**2-qs(2)**2))

	fsnext(2) = -(qs(2)+lambda*(qs(2)**3+qs(2)*qs(1)**2)- &
                2*lambda1*qs(1)*qs(2))
 
             
		pupdate(1) = 0.5*(fs(1)+fsnext(1))
		pupdate(2) = 0.5*(fs(2)+fsnext(2))
		ps(1) = ps(1) + pupdate(1)*dt
		ps(2) = ps(2) + pupdate(2)*dt
		fs = fsnext
    enddo
	ef = (qs(1)**2 + qs(2)**2 + ps(1)**2 + ps(2)**2)/2. + &
              0.25*lambda*(qs(1)**2 + qs(2)**2)**2 + &
              1./3.*lambda1*(qs(1)**3-3*qs(1)*qs(2)**2)

	momf = qs(1)*ps(2)-qs(2)*ps(1)
!	check energy conservation
!	print *, "change of energy during MD:", ef-e0
!	print *, "change of momentum during MD:", momf-mom0
    return
  end subroutine mdvv

end module md_vv
